Collisions of anisotropic two-dimensional bright solitons in dipolar Bose-Einstein 

condensates 
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We investigate the coherent collision of anisotropic quasi-two-dimensional bright solitons in dipolar 
Bose-Einstein condensates. Our analysis is based on the extended Gross-Pitaevskii equation, and we 
use the split-operator method for the grid calculations and the time-dependent variational principle 
with an ansatz of coupled Gaussian functions to calculate the time evolution of the ground state. We 
compare the results of both approaches for collisions where initially the solitons are in the repelling 
side- by-side configuration and move towards each other with a specific momentum. We change the 
relative phases of the condensates, and introduce a total angular momentum by shifting the solitons 
in opposite direction along the polarization axis. Our calculations show that collisions result in 
breathing-mode- like excitations of the solitons. 

PACS numbers: 03.75.-b, 05.45.-a, 67.85.-d, 34.50.-s 



I. INTRODUCTION 

Bose-Einstein condensates (BECs) of magnetic atoms 
have attracted much attention since their experimental 
realization with ^^Cr atoms [1 . Recently, the creation of 
condensates of ^^^Dy [2] and ^^^Er [3| atoms with much 
larger magnetic moments than ^^Cr have also been re- 
ported. Furthermore, there has been fast progress to- 
wards the condensation of polar molecules with electric 
dipole moments [4 , where the dipole-dipole interaction 
(DDI) is even more dominant. A review of the physics of 
dipolar bosonic quantum gases has recently been given by 
Lahaye et al. [5 . The features of the DDI being a non- 
local long-ranged and anisotropic interaction give rise to 
a variety of new effects. One example is the creation of 
solitary waves, where in analogy to nonlinear optics the 
effects of dispersion and nonlinearity may cancel each 
other. This leads to a condensate with a shape constant 
in time. The experimental realization of one-dimensional 
solitons in self-attractive BECs of ^Li atoms has been 
reported [6^. Tikhonenkov et al. have theoretically pre- 
dicted 2D solitons [7] and Koberle et al. have proposed a 
realistic experimental setup for the creation of a 2D soli- 
ton [HI. An exciting aspect of multidimensional solitons 
is their anisotropic nature, based on the in-plane polar- 
ization of the dipoles of such solitons. 2D solitons have 
already been studied using a variational ansatz with a 
single Gaussian and with coupled Gaussian functions [9] . 
Adhikari et al. have recently investigated axially sym- 
metric and vortex solitons on a one-dimensional optical 
lattice [1^. Note that in contrast to systems with har- 
monic traps, where the density distribution in the trap 
direction is an approximate Gaussian, systems in an op- 
tical lattice will have an exponential density distribution. 



The collision of axially symmetric bright 2D solitons 
has been studied by Pedri et al. [Ij and Adhikari et 
al. [10 . Pedri et al. investigated a system with dipoles 
aligned parallel to the harmonic trap, while Adhikari et 
al. used an optical lattice instead. In both cases, the 
sign of the DDI has to be inverted by fast rotation of the 
orientation of the dipoles [12 . The resulting interaction 
energy becomes Ud{R) = — a(3cos^i? — l)/i^^, where 
is the angle between the polarization axis and R = r — r' . 
The factor a can continuously be changed from —1/2 to 
1. This provides the possibility to change the dipolar 
interaction from attractive to repulsive. 



In addition to the analysis of the collision of 2D soli- 
tons. Young et al. [13] have investigated the collision of 
one-dimensional bright and vortex solitons. The inves- 
tigations in [THS] concentrated on the creation and the 
stability of 2D solitons with respect to small perturba- 
tions. However, one important property of solitons is 
that their shape is constant in time even when they are 
moving. Therefore, the collision of two solitons is an ad- 
equate scenario for the investigation of soliton dynamics 
far beyond small excitations. The influence of the non- 
linear contact interaction and the DDI are of particular 
interest in such calculations. 



As mentioned above, the creation of a BEC of magnetic 
atoms has been realized with a variety of species. Our 
results are valid for all dipolar systems, but we will add 
the corresponding values for a system with 20 000 ^^Cr- 
atoms per soliton in parentheses. 



At sufficiently low temperatures, the dynamics of a 
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H{t)'^{r, t) = (-A + Fhar + + Vi) *(r, t) 

=idt^ir,t), (1) 
with T4ar =7'y' , = Sira |*(r,t)|' , 



V'(r-,i + At) =(r|C/(At)|V') 



r — r'\ 



Here a is the scattering length and ^ designates the 
mean-field wave function. The dipoles are polarized along 
the z-axis, so that ^ is the angle between the z-axis and 
the vector r — r'. We choose the ^-direction as the axis of 
confinement perpendicular to the polarization axis where 
= 20 000 (420 Hz), while the condensate is free in x- 
and 2:-direction. All simulations deal with condensates of 
low densities, and only a small period of time in which the 
two condensates merge to one transient condensate with 
higher density. This means that we do not need to take a 
three-body-loss term [8 into account, as the resulting ab- 
sorption images integrated along the ?/-axis) would 
only be slightly affected. We checked this assumption for 
the calculation of the collision without difference in phase 
and without angular momentum which up to the time of 
t = 0.06 {t = 0.001 corresponds to 15 ms) resulted only 
in a loss of about 5.5% of the particles. 

As has been shown in [9 , solitons only exist in a cer- 
tain range of values of the scattering length, which can 
be tuned by the use of Feshbach resonances [15 . For too 
large values, the condensate will disperse, while too small 
values lead to the collapse of the condensate. In the fol- 
lowing the scattering length is chosen to be 0.14 (12.7aB, 
where ae is the Bohr radius) . 



II. NUMERICAL APPROACH 



The main theoretical task for the grid calculations is 
how to apply the time evolution operator U = e~*^^ on 
a state For this, one splits U symmetrically by using 
the Baker-Campell-Hausdorff formula |16J 



— 7-TAf - 



iVAt -i^TAt 



(2) 



where V = Vhar + V'sc + Vd- One projects the action of 
the approximated time evolution operator on the basis of 
the position operator and makes use of the possibility to 
insert / diy\u) = 1: 
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^-ip'r' ^-iVir')At^ir'p^-i4At^l^j^y ^3) 

The structure of ([3| suggests the following algorithmic 
procedure: 

• Fourier transform of '^(r) in order to obtain ^p{p) 

• Multiply by e'^^^i 

• Inverse Fourier transform to real space 

• Multiply by e-^^^^)^* 

• Fourier transform to momentum space 

• Multiply by 6"^^ 

• Inverse Fourier transform to real space 

The potential V consists of the harmonic potential, the 
scattering potential and the DDI potential. The scatter- 
ing potential and the DDI potential have to be calculated 
at each time step. The latter can be evaluated by means 
of the convolution theorem, which results in two more 
Fourier transforms: 



(4) 



Here k and kz denote the momentum and the momen- 
tum in z-direction, respectively. Altogether, we have to 
perform six Fourier transforms for each time step. Note 
that the first and last Fourier transforms described in 
the algorithmic procedure of the time evolution are only 
necessary if one is interested in physical quantities whose 
evaluation requires the wave function in real space. 

For the simulations, the spatial domain was discretized 
with up to 512 X 128 x 512 grid points. Since this scheme 
is numerically very demanding, it has been implemented 
for graphics processing units (GPUs) using CUDA, en- 
abling a very high degree of parallelization. Using the 
Tesla C2070 improves the performance of our algorithm 
by a factor of about 80 for double precision in compar- 
ison to the corresponding C algorithm using the well 
known FFTW library for computing the discrete Fourier 
transform on a IBM System x3400 with a Quad-Core In- 
tel Xeon Processor E5430 (2.66GHz 12MB L2 1333MHz 
80w) and 4 x 4GB PC2-5300 CL5 EGG DDR2 Chipkih 
Low Power FBDIMM 667MHz. 



3 



To investigate the coherent cohision of sohtons we have 
apphed the fohowing procedure. The first step is the 
computation of the ground state of one condensate using 
the spht-operator method with imaginary time evolution 
{t = —ir). Afterwards we double the size of the grid 
in the x-direction and place two solitons in the repelling 
side-by-side configuration. The distance between the con- 
densates is chosen such that they do not feel the mutual 
dipole-dipole interaction. To introduce momentum in the 
system, we multiply the left hand-side of the wave func- 
tion by a plane wave e*^^ (for the soliton moving to the 
right) and the right hand-side by e~*^^ (for the soliton 
moving to the left), respectively. 



III. TIME-DEPENDENT VARIATIONAL 
ANSATZ 

Variational calculations using coupled Gaussian wave 
packets (GWPs) have shown to be a full-fledged alterna- 
tive to numerical grid calculations for the calculation of 
ground states of dipolar BECs [9,T7l. The applicability 
of such ansatzes to dynamical simulations is a challenging 
task. The decisive extension of the previous work [9| dZ] 
is that additional translational and rotational degrees of 
freedom are included in the ansatz with coupled GWPs to 
describe the dynamics of the condensate wave function. 
For the convenience of the reader we shortly review the 
time-dependent variational principle (TDVP) in this sec- 
tion, and subsequently apply it to the ansatz of coupled 
GWPs. We make use of the TDVP in the formulation of 
McLachlan flS^ where 6 is varied such that 



/= ||i^-i^^(t)||^ =min. 



(5) 



and set = ^ afterwards. The wave function ^ is con- 
sidered to be parametrized by the variational parameters 
^ = ^(z(t)). The minimization of the quantity / in 
Eq. ^ leads to 



which can be written in the short form 
Kz = -ih, 



(6) 



(7) 



with the positive definite Hermit ian matrix K. We use a 
linear superposition of Gaussian wave packets (GWPs) 

k=l 

-E/, (8) 



k=l 



as an ansatz for the wave function in Eq. ([s]). In gen- 
eral, are 3x3 complex matrices (determining the 



width and the orientation of the GWP), and are 
three-dimensional real vectors (representing momentum 
and center of the GWP) and 7^ are complex numbers 
(where the real part stands for the amplitude and the 
imaginary part for the phase of the GWP, respectively). 
In this work we will make use of the strong confinement 
in one direction perpendicular to the dipole axis and omit 
the translational and rotational degrees of freedom in y- 
direction 



0, pS = o, 



0, 



(9) 



with a = x^z. Inserting the ansatz Eq. ^ in Eq. (|6]), 
sorting the result by powers of x and identifying these 
terms with the coefficients of a time- dependent effective 
harmonic potential 



eff 



• V-{X ■ 



XV2 X , 



(10) 



yields the equations of motion (EOM) for the variational 
parameters 



-4i {A'' 



P 



(11a) 

'(lib) 
(11c) 

(lid) 

If we write Eq. ^ explicitly for GWPs, the set of linear 
equations for z can be rewritten to one for the vector v 
containing the coefficients of V^^ 



= -Rev^-2 Im A^ {q^ - 2p^) - 2 Re V^^q^ 

= + i {ReA^y^ {Imv^ -^2lmV2^q^) 
= 2iTr - iq^Vi'q^ + Ap^A^q^ + i {p^Y 



iq^p^ 



Kv = r , 



(12) 



for details, see p!9l[20|. With the transformation given 
in Appendix |A] the EOM can now be integrated with a 
standard algorithm such as Runge-Kutta, where Eq. (12) 
has to be solved at every time step. The right-hand side 
vector r with the components 



N 



r' = Y.{9'\x^x-pV{x) \g'^) 



(13) 



where / = 1,...,A/'; = 1,...,3; < n -\- m < 2, 

contains integrals of the potentials in the GPE. It is one 
of the most important advantages of the method that 
nearly all of these integrals can be calculated analyti- 
cally. However, the dipolar integral (^| |^) can only 
be calculated analytically for GWPs centered in the ori- 
gin without the additional translational degrees of free- 
dom introduced in the ansatz ([8|. The analytical and 
numerical treatment of the dipolar integral is presented 
in Appendix [B) 

The procedure for the calculations is as follows: At 
first the equations of motion (11) for one soliton are in- 



tegrated in imaginary time, with the wave function be- 
ing normalized after every time step. Afterwards every 
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GWP of the wave function is copied and the resulting 
two sohtons are positioned in the same way as given in 
Sec. [TTj Then for each GWP a corresponding momentum 
= ±p^ex is added. For this starting configuration the 
EOM are finahy integrated in real time. 



IV. RESULTS 

In Fig. [T] three grid calculations of colliding solitons 
without angular momenta prepared in the way given 
above are shown. For no phase difference constructive 
interference occurs and the condensates merge and split 
up in two solitons again. Note that the condensates after 
the split-up {t = 0.049, t = 0.001 corresponds to 15 ms) 
have a larger spatial distribution than before {t = 0.011). 
This indicates that the transfer of kinetic energy to in- 
ternal energy has excited the solitons. This might either 
induce the dispersal of the solitons or lead to breathing- 
mode-like oscillations. The column in the middle shows a 
simulation with a difference of ^ = 7r/2 in phase, result- 
ing in a collision where the soliton on the right eventually 
has a lower amplitude than the one on the left, so that we 
do not have symmetric behavior anymore. The transfer 
of kinetic energy is not as large as for = 0, resulting in 
only slightly larger condensates at t = 0.049. In the case 
of a collision with a difference of = tt in phase we can see 
destructive interference (column on the right), the soli- 
tons effectively repel each other. The transfer of kinetic 
energy into internal energy is once again smaller, corre- 
sponding to a just slightly larger condensate at t = 0.049. 
The occurrence of the broken symmetry in x-direction 
can be understood if one considers that a difference in 
phase of (/) = and (j) = tt yields a wave function which 
is an eigenfunction of the parity operator, in the sense 
of II±^(r,t) = ±^(-r,t). A difference of ^ = 7r/2 on 
the other hand does not result in an eigenfunction of the 
parity operator, thus yielding an asymmetric dynamic of 
the condensates. 

In Fig. [2] we compare the results for grid calcula- 
tions and the variational ansatz for simulations, where 
we shifted both condensates in opposite directions along 
the polarization axis in order to introduce angular mo- 
mentum. Both approaches are in very good agreement 
with only slight differences, in particular for times where 
both condensates merge, and when comparing the ex- 
tensions of the solitons at t = 0.049. It is remarkable 
that a total number of only six GWPs is sufficient to re- 
produce the structures of the grid calculations and give 
the correct result for the configuration at the end of all 
three simulations. The first case without a difference in 
phase (Fig.|2^) once again leads the solitons to merge and 
split up afterwards, while a transient eddy-like structure 
appears in the course of the collision. As in the case 
with no angular momentum, the solitons either seem to 
disperse, or a breathing-mode-like oscillation has been 
excited. The amount of kinetic energy which has been 
transferred is lower than in the former case, leading to 
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FIG. 1. (Color online) Absorption images (IV^I^ integrated 
along the ^-axis) of grid calculations for the collision of soli- 
tons. The value of the momentum for each soliton is A: = 10 
(velocity v = 127/im/s) and the field of view is 1.4 x 1.4 
(135/im X 135/im). All absorption images have been normal- 
ized to the maximum value. Left column: Absorption images 
for a collision without difference in phase. Middle column: 
Absorption images for a collision with a difference of = 7r/2 
in phase. Right column: Absorption images for a collision 
with a difference of = tt in phase. 

condensates with smaller extension at t = 0.049 than 
their corresponding condensates in the simulation pre- 
sented above. A difference of (/) = 7r/2 in phase (Fig. [2]3) 
shows a similar behavior as in the case without angular 
momentum, resulting once again in an asymmetric sit- 
uation where after the collision the condensates do not 
have the same amplitudes anymore. But for finite an- 
gular momentum one may actually speak of a merged 
condensate at t = 0.031. Finally the collision with a 
difference of = tt (Fig. ^) shows the condensates effec- 
tively repelling each other, but in this case introducing 
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(a) (b) (c) 




FIG. 2. (Color online) Absorption images for grid calculations and the variational ansatz of the simulation of two colliding 
solitons with angular momentum. All absorption images have been normalized to the maximum value. The parameters are the 
same as given in Fig.[l] The columns (a), (b) and (c) show calculations for a difference of = 0, (j) — n/2 and = tt in phase, 
where the left column is the result of the grid calculations and the column on the right hand-side presents the results of the 
variational ansatz. For all three calculations six GWPs (three for each soliton) were used. The variational calculation is able to 
reproduce the transient ring-like structure during the collision for a difference of = tt in phase and yields the correct results 
for the configuration at the end of all three calculations. 



angular momentum leads to a transient ring-like struc- 
ture. The extension of the condensates after the collision 
is much larger compared to the case with no difference 
in phase, which means that the amount of transferred 
kinetic energy in internal energy is larger than in the for- 
mer case. The case with angular momentum is suited 
best to show how the transfer of kinetic energy affects the 
spatial distribution of the condensate. In Fig. [3] we show 
the kinetic energy as a function of time for the collisions 
with angular momentum. Comparing the curves in Fig. |3] 
with the absorption images in Fig. |2j it is obvious that a 
larger transfer of kinetic energy implies a larger conden- 
sate at t = 0.049. The slightly smaller transfer observed 



at the end of the full-numerical calculations (this leads 
to a larger extension of the solitons after the collision c.f. 
Fig. [2| originates from finite grid sizes and thus has no 
physical meaning. Variational calculations show an os- 
cillation of the kinetic energy for large timescales, which 
corresponds to the excitation of the solitons. 

The amount of kinetic energy transferred into inter- 
nal energy of the solitons depends on the overlap of the 
wave functions during the collision process. A large over- 
lap of the solitons enhances the nonlinear coupling in the 
GPE as |^(x,t)p increases and a small one diminishes 
the coupling. This can be seen best in Fig. [l] (right col- 
umn) where the destructive interference for the calcula- 
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FIG. 3. (Color online) Kinetic energy as a function of time 
for the collisions with angular momentum (Fig. [5]). The dots 
show the numerical results, the lines show the results obtained 
by the variational calculations. The kinetic energy increases 
while the condensates merge. After the split up, the conden- 
sates have a lower kinetic energy than before, indicative of a 
transfer from kinetic to internal energy, thus resulting in exci- 
tation of the condensates. The inset shows the kinetic energy 
obtained by the variational calculations for large timescales. 
The oscillatory behavior indicates the excitation of the soli- 
tons. Note the larger kinetic energy obtained by the grid 
calculations shortly after the split up. This is due to the fi- 
nite grid size, which manifests itself in oscillations of the wave 
function's amplitude for large times. 



tion with phase difference ^ = tt leads to |^(0,t)p = 0. 
For the corresponding calculation with nonzero angular 
momentum (Fig. [2]3) we find |^(0,t)p = 0, too. How- 
ever, the ring-like structure increases the overlap during 
the collision. 

In Fig. [i] the variance A^r = (cr^) — (a)^ with cj = x^z 
is plotted as a function of time. The variance has been 
calculated for the three GWPs representing the solitons 
on the left-hand side in the starting configuration and for 
the GWP which has the largest amplitude after the col- 
lision process. This dominant GWP shows oscillatory 
behavior while the other GWPs with much smaller am- 
plitudes describe particles leaving the soliton. This effect 
can hardly be seen in the absorption images in the upper 
panel of Fig. [4j However, the absorption images show 
that a soliton still exists, although this would be difficult 
to see in an actual experiment due to the very long time 
scale. 

We have also performed simulations with smaller and 
larger momenta of the solitons. The former case leads 
to one merged condensate which does not split up again 
after the collision but shows oscillatory behavior. This 
is very similar to the collision presented in [11 . In the 
latter case the wavelength of the interference pattern is 
smaller and becomes more pronounced. Note that grid 
calculations with high momenta are problematic, because 
the condensates quickly reach the edge of the grid. An 
approach with a variational ansatz is better suited to 
analyze these scenarios. 
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FIG. 4. (Color online) Variance A^j of a single soliton (right 
soliton in the absorption images in the upper panel) as a func- 
tion of time. The collision occurs at t ~ 0.03. The thick solid 
line and the double-dashed line show the A^ and A^ variance 
of three GWPs, respectively. The thin solid line and the nor- 
mal dashed line show the A^ and A^ variance of the dominant 
GWP g^. 



V. CONCLUSION 

We have studied the collisions of anisotropic two- 
dimensional bright solitons in dipolar Bose-Einstein con- 
densates both with a fully-numerical ansatz and a time- 
dependent variational principle with coupled Gaussians. 
The calculations presented show that the collision pro- 
cess leads to an energy transfer from kinetic energy to 
"inner" energy of the solitons which leads to excited soli- 
tons with larger extent. The absorption images show very 
good qualitative agreement of the results gained by the 
two different methods. 

The advantages of the grid calculations are the sim- 
plicity of the numerical scheme (although the implemen- 
tation for the massively parallel computation requires 
some effort), the freedom in describing all different shapes 
of wave functions, and the numerical stability of the 
method. The advantages of the variational calculations 
are the much smaller numerical effort, enabling one to run 
long calculations on standard PCs, the independence of 
finite grid size, and the small amount of parameters to 
be saved. 

Both methods can be used to simulate the time- 
dependent GPE, supporting each other mutually. One 
further application would be the inclusion of additional 
external potentials such as optical lattices and the com- 
parison of the methods in such scenarios. Our results 
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should stimulate experimental efforts to study the colli- 
sions of 2D anisotropic solitons. 



VI. ACKNOWLEDGEMENTS 

We thank Boris Malomed for valuable discussions. 
This work was supported by Deutsche Forschungsgemein- 
schaft. R.E. is grateful for support from the Landes- 
graduiertenforderung of the Land Baden- Wiirttemberg. 



Appendix A: Transformation to CB-variables 



The direct numerical integration of Eq. (11) leads to 



numerical difficulties [21 . These can be dealt with by the 
introduction of two auxiliary matrices B and C. With 
A = BC~^ the equations of motion for the width matri- 
ces can be written as 



-4i {A^y ^iV^^, 



(Al) 
(A2) 



where C and B are 3x3 complex matrices. Omitting 
the index k we obtain from these equations 

B-^BC-^ - C-'^C = -AiB-^A^ + 15" V2 (A3) 
^ B-^B -C-'^CC = -4iC-^B ^iB-W2C . (A4) 

By comparison we yield the equations of motion for C 
and B 



B^ = iV^C^ , 
^k ^ ^ 



(A5) 
(A6) 



The reduction ^ can be done for those matrices, too. 
Note however, that the matrices B and C do not preserve 
the same symmetry as the matrices A which are complex 
symmetric. Therefore, all five complex entries in B and 
C have to be integrated. 



Appendix B: Solution of the dipolar integral 

The calculation of the dipolar integrals needed in the 
TDVP {^\ a^/3^1/d 1^) with a, /3 = x, ^, z and < n + 
m < 2 is shown here for the simplest case n = m = 0. 
The other integrals are calculated analogously. We start 



from the six-dimensional non-local integral 
{<S\Va\^)= jj dhdh' g'\r)gi\r')g\r')g\r) 

l,k,j,i 

^ 3{z-z')^\ 1 



X 1 



\r — r'f I \r — r'f 



■ (Bl) 



By the use of the convolution theorem of Fourier analysis 
we can evaluate one of the three-dimensional integrals 
directly, while the inverse Fourier transform 



(B2) 



remains to be done. Here /q^ 



denotes the overlap integral 
of the Gaussian functions k and /, and we have used the 
abbreviations 



j^klij 
pklij 



-Ipkl 



(B3) 
(B4) 



with = + and = p_f and analo- 
gously for i and j. The integral ( |B2[ ) can be split in 
two parts, one leading to a shift in the scattering length 
(this is the short-range part of the DDI) and a second 

jkl pj jklij 
l^Lk.jA ^0 ^0 ^2 



I^^ J^'^^^ . After a principal 



part (^|yd,eff 1^) 
component analysis of the exponential in Eq. (|B2[) the 



analytical integration in /c^-direction is possible when we 
make use of Eq. ^ . The remaining result reads 



00 



with Xc — ^/(\^rx)J2^ Xs = Y^(l — x)/2, the coeffi- 
cients Co, ci of the rotation matrix from the principal 

component analysis and the Faddeeva function w{z) = 
2 

e~^ erfc(— iz). The numerical evaluation of this integral 
can efficiently be performed by a Taylor expansion of the 
Faddeeva function for which the single terms can be ob- 
tained by a recursion formula and using a Chebyshev 
quadrature for the x-integration. To improve the result 
we apply a Fade- approximation to the Taylor series. 

The numerical integration of the dipolar integrals is the 
crucial part in this method. Dependent on the number of 
Gaussian functions there is a total number of Cnum = 
{N^ + N'^ — 2N)/A integrals to be calculated numerically 
and Ceiiiptic = N{N + l)/2 which can be expressed in 
terms of elliptic integrals. 
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